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ABSTRACT 



The classical problem of spherical homologous gravitational collapse with a 
polytropic equation of state for pressure is examined in Lagrangian fluid coordi- 
nate, where the position of each initial fluid element r\ — r(0) is followed in time 
by the evolution function y(t). In this Lagrangian description, the fluid veloc- 
ity v = dr/dt = rjdy/dt is not a fluid variable, contrary to the commonly used 
Eulerian fluid description. As a result, the parameter space is one dimensional 
in 7], in contrast to the (x, v) two-parameter space of Eulerian formulation. In 
terms of Lagrangian coordinate, the evolution function y(t), which is not lim- 
ited to a li near time sca ling, agrees with the well established parametric form 
of Mestel ( jMestel 1 119651 ) for cold cloud collapse. The spatial structure is de- 
scribed by an equation which corr esponds to the one derived by Goldreich and 
Weber (IGoldreich &: Weber I Il980l ). The continuous self-similar density distri- 
bution presents a peaked central core followed by oscillations with decreasing 
ampl itude , som ewhat reminiscent to the expansion-wave inside-out collapse of 
Shu (jShu Ill977l ). This continuous solution could account for the planetary sys- 
tem of a protostar. There is also a disconnected density distribution, which could 
be relevant to cavity formation between the highly peaked central core and the 
external infalling envelope of a magnet ar- in- a- cavity pre-supernova configuration. 



Subject headings: gravitational collapse 
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1. Gravitational Collapse 



The fundamental issue of gravitational collaps e of a large but finite massive gas cloud 



concerns the star and nebula formation. Bonner (IBonner 



19561) had demonstrated that 



the equilibrium of a finite i sothermal gas cloud under its self-gravity could be unstable as 
the mass increases. Mestel (IMestel Ill965l ) solved the collapse equation of a cold cloud with 
v = dr/dt as the trajectory velocity, which is the Lagra ngian fluid repre sentation, to ob- 
tain the celebrated self-similar parametric solution. Lin ( ILin et al. Ill965l ) showed that an 
oblate cold spheroid would evolve to a disk and a prolate o ne would evolve to a spindle. 



Bodenheimer and Sweigart ( iBodenheimer fc Sweigart I Il968l ) studied numerically the evo- 
lution sequence of a collapsing gas cloud with finite preesure again with v = dr/dt under 
different i nitial density distributions and different surface boundary conditions. Penston 
19691 ) analyzed analytically the Lagrangian cold collapse with a smooth maximum 



flPenston 



for density, and also put forward a self-similar analysis of an isothermal collapsin g sphere 



in Eu lerian description, where the fluid velocity v as one of the variables. Larson (ILarson 



19691 ) examined numerically through a set of conservation equations of mass, momentum, 
and energy in Eulerian representation the format ion of a p rotostar, and presented in ap- 



pendix C the isothermal similarity solutions. Shu (IShu 



19771 ) studied anew the homologous 



collapse of an isothermal sphere with the Eulerian conservation equations. He interpreted 
the singular solutions, where the coefficients of the nonlinear differential eq uations vanish , 
and constructed the expansion- wave, inside-out collapse scenario. Hunter (iHunter 1 119771 ) 
added a ne w class of isothermal self-si milar solutions on previously known ones. Goldreich 
and Weber (IGoldreich fc Weber Ill980l ) examined the homologous collapse of a stellar core 
with a polytropic equation of state, and arrived at a differential equation that describes the 
radial structure. Perturbation anal ysis was used to study the nonsp herical modes of the stel- 
lar core. Whitworth and Summers (IWhitworth fc Summers Ill985l ) brought to the attention 
the importance of the stability of the initial isothermal gas cloud and the external driving 
pressure. These two factors transform each known solution into a continum. 



Recently, there h aye been renewed effo rts on gravitati onal collapse. Fatu zzo (iFatuzzo et al 



20041 ) , Lou and Wang (ILou fc Wang 1120061 ) , Lou and Gao (ILou k, Gao 1120061 ) have considered 
self-similar solutions in Eulerian form with polytropic equation of state, and pr esented inno 



yated solutions for astrophysical phenomena. Furthermore, Lou and Shen (jLou &: Shen 



20041 ) put forward the envelope-expansion (EE) and core-collapse (CC) solution for an 
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isothermal cloud. This EECC solution is to mo del outflowing stellar winds while the red 
giant collapses. Bian and Lou ( Bian fc Lou 1 2005 ) examined the shock flows of an isothermal 
sphere to better understand the shock structure. 

We remark that, except cold gas clouds, clouds with isothermal or polytropic equation 
of state have been treated under Eulerian fluid description where velocity v = v(r,t) is 
one of the variables. The time dependence of each variable is constructed as some specific 
power of time t. Th e spatial dependence is given by a set of nonlinear differential equations 
( iLou fc Shen II2004J ). The singular solutions of this system define the singular sonic surfaces, 
which separate subsonic and supersonic regions in the (x, v) two-parameter space. Along the 
sonic surfaces, the flow is subsonic/supersonic relative to the similarity x profile. Collapse 
solutions cross from subsonic to supersonic regions along the sonic lines. Sometimes shock 
waves are needed to bridge the crossing to meet the boundary conditions. Aimed to describe 
the collapse, representing the time dependence with different powers of a linear time scaling 
for different variables would certainly fail to describe the time history of collapse. We could 
only expect such self-similar solutions with a linear time scaling be valid during a specific time 
span and over a specific radial interval. Here, we use the alternative Lagrangian description, 
where the fluid velocity is not a variable, to construct the self-similar collapse having a 
polytropic equation of state. The evolution function is solved consistently together with the 
spatial structures. 
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2. Lagrangian Similarity 



Gravitational collapse of a spherical cloud is described by 



§ + v-M = , (1) 
'^V(.-V)^ = 0. (3) 



dt VP 7 / VP 7 . 
Here, p is the mass density, v is the gas cloud velocity, p is the gas pressure, 7 is the polytropic 
index, M* is a point mass at the center, if any, and M(r, t) is the gas mass within a sphere 
of radius r at time t where 

M(r,t) = [ 4vrr /2 p(r' ' ,t)dr' . (4) 



Jo 

We note that Eq. Q] can be written alternatively as 

DM dM 

_ + l ,_ = C W = 0, (5) 

where C(t) is a constant independent of r and varies only in time. By taking C(t) = 0, 
we would have an inhomogeneous cloud that is stationary in time. For a spherical collapse 
with v = vf, we use Lagrangian fluid representation with r(t) = r(0)y(t) = i]y(t), where the 
Lagrangian label 77 is the initial position of the fluid element, and y(t) is the dimensionless 
evolution function with y(0) = 1 that describes the time history of the fluid element. The 
velocity can then be written as 

d y (r\ 
v = rj- . (6) 

Contrary to the Eulerian fluid description, the fluid velocity v here is not a variable which 
substantially simplifies the analysis. 

We now transform the Eulerian independent variables (r(t),t) to Lagrangian indepen- 
dent variables (rj,y(t)). Furthermore, we note that the total time derivative in Eulerian 
representation (d/dt + vd/dr) corresponds to the partial time derivative in Lagrangian rep- 
resentation d/dt. We now write the variables in saparable form of y and 77, and try to 
determine the evolution fun ction y(t) and the spa tial configuration in 77 under similarity. 



Following Tsui and Serbeto (ITsui fc Serbeto 



20071 ) , we have from Eq. \5\ and Eq. [3] 



M(r), y) = M ^M( V ) , (7) 
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F — — — ^-F(n) (8) 

where M carries the physical dimension of M(rj, y) such that M{rj) is dimensionless, and 
likewise are Po/Po an d F(t]). We note that F is only a function of entropy which remains 
constant in time and in space. For this reason, we have F(rj) = 1. Besides, from Eq. [Hand 
with p(r),y) = p Q Y(y)p(r}), we have 

M(r,t) = y 3 f ^T/ 2 p Q Y{y)p{rf)dTf . (9) 
Jo 

Comparing with Eq. [7| gives 

p(v,y) = po-^p(v) , (10) 

y.i 

M(z) = [ 3z' 2 p(z')dz' , (11) 



Mo = ^ , (12) 

where we have introduced the scaling a to write 77 in normalized form z = at]. Making use 
of Eq. [SJ we now have 

p(v,y) = Po-kzfi.n) • (13) 

As for Eq. [2j in terms of Lagrangian coordinate and with 7 = 4/3, it reads 

2 d 2 y Po a 2 lldp GM*a 3 GM a 3 - 

V W = —^-pld-z ~ 3 — M = C - (14) 

This equation has the temporal and spatial parts separated with C as the separation constant. 
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Evolution Function 



Choosing C = —NGM a 3 in Eq. [[Hand denning r = (GM a 3 ) 1 / 2 t = cot, the temporal 
part of Eq. [TH is 



fy 

dr 2 



N 



(15) 



Taking N = would give dy/dr = —n and y{r) = y(0) — nr, which reproduces the linear 
time scaling of the Eulerian formulation. For n = 0, we have the equilibrium state of a 
gravitating shherical cloud. For N > 0, y(r) has a downward curvature. This is compatible 
to an inflowing fluid to the center with y{r) decreasing from its initial value of y(0) = 1. 
The first integral of this equation is 

-H . (16) 



dy 
dr 



2N\ 

y J 



The terms on the left side could be interpreted as the kinetic and potential energies and the 
constant —H on the right side could be regarded as the total energy. The inflow rate is given 
by 



dy 

d:T 







1/2 




± 


Kt-*): 




Kt-")] 











1/2 



(17) 



With H = 2N, the initial velocity dy/dt = starts from zero, and Fig.l shows the evolution 
function y{r) and its time derivative dy/dr respectively. 

We remark th at this evolut ion function actually corresponds to the parametric time 
function of Mestel ( Mestel Hl965l l. Taking y = cos 2 9, Eq. [[6] with 2N = H gives 



dy 
dr 



Htan 2 6 



With dy/dr = 2 cos 9 sin 9d9 / dr, we have 

2 cos 2 9d9 = (l + cos29)d9 = H 1/2 dr 
Integrating once then gives 

9 + - sin 20 = H 1/2 r , 
2 



18) 



which is the well established parametric solution of Mestel. 
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With C 



4. Spatial Structure 

—NGM a 3 , the spatial part of Eq. [Til is 
^ 2 1 11 dp M 



---^ + -7 = N 



z* 



Po GM a 3 p z dz 

With p = p 7 , C 2 = jpo/po, and vjf = 2GM Q a where // denotes free fall, we get 



CI 



J ff 



(7 - 1) dz 



2 9 -7-1 

z —p 1 



p(z')3z' 2 dz' + Nz 3 



Denoting a = C 2 /vjf, we get 



« 2 9 
2aZ d~z P 



7-1 



-3(7-1) / (p(z') - N)z l2 dz' 



(19) 



(20) 



(21) 



(22) 



With 7 = 4/3, and writing p 1//3 = q and p = q 3 , Eq. [2T1 becomes 

2az 2 ^- = - f \q 3 {z')-N)z' 2 dz' . 
dz Jo 

Differentiating this equa tion once leads to the follo wing; equation, which is similar to Eq.6 
of Goldreich and Weber (IGoldreich fc Weber Ill980l ). 



2a- 



1 d 



,dq 



+ Q 



N 



(23) 



z 2 dz \ dz 

In Goldreich and Weber, iV is their free parameter A which approaches zero, and g(0) is chosen 
as unity. In our case, a and q(0) are the fr ee parameters. We remark that Bodenheimer and 
Sweigart (jBodenheimer &: Sweigart Ill968l ) had pointed out the importance of the parameter 
a. The most important difference is our constant N — 1, which comes from the separation 
constant C. Should we take C = and thus iV = 0, we would have dy/dr = —n and 
U = 2/(0) — nT from Eq. [JSJ This would recover the linear time scaling of the Eulerian 
similarity and Eq. [23] would equal to Eq.6 of Goldreich and Weber. The fact that N = 1 not 
only provides an evolution function of Mestel, but also changes qualitatively the nature of 
spatial solutions. 

With a = for a cold fluid, the equilibrium requires dy/dr = with y(r) = y(0) = 1. 
This means N = and n = 0, and it is obvious from Eq. [22] that such a cold fluid equilibrium 
does not exist, other than q(z) = 0. As for the collapse solution of iV > 0, we have the 
uniform density sphere solution of 



p{z) = q 3 {z) = N 
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This states that a uniform density cold sphere under homologous infall stays as a uniform 
sphere. By taking N = 1, po amounts to the density of the sphere. Making use of Eq. [TTj 
M(z max ) = z 3 max , and therefore Mqz^^ is the mass of the entire uniform sphere of z max . 

With a / 0, equilibrium configuration with N = gives a monotonically decreasing 
profile that ends at a finite z max < oo with q(z max ) = 0. With q(Q) = 2 and a = 100, the 
profiles of p = q 3 and q are shown in Fig. 2 with z max = 40. As for the collapse solution, we 
note that when q > N — 1, the integral of Eq. [22] is positive, and dq/dz is negative with 
g decreasing. As the upper limit z of the integral increases, q would decrease to less than 
unity, and the integral would begin to decrease making dq/dz to turn around. We would 
then get oscillating solutions about q — 1, as shown in Fig. 3. This oscillating structure with a 
progressively smaller amplitude shows that the central core collapses faster than the external 
oscillating envelope. This spatial str ucture bear s great ressemblance to the expansion-wave 
inside-out collapse scenario of Shu (IShu I Il977| ). For a finite system with z max < oo, the 
solution is a function of q(0) for a given a. As q(0) increases with the same a, the first 
minimum becomes more extensive, and followed by damped ripples about q = N = 1. The 
peak about q(0) corresponds to the central core, and the oscillating part amounts to the 
external envelope. 

This spherical oscillating envelope structure of q(z) could be very relevant to protostar 
planetary system. The locations of the gaseous planets would be given by the oscillating 
peaks of p = q 3 . In the presence of rotation, as collapse proceeds, the spherical cloud would 
flatten to the equatorial plane, and the high density shells in the spherical cloud would 
become dense rings in the equatorial disk forming gaseous planets. Rotation would also 
stabilize the planets on their orbits as the central pre-protostar continues to collapse. Since 
the peaks are decreasing in amplitude, the gaseous planets would have a decreasing mass 
as distance increases. This happens to be the case of the giant planets of our Solar system, 
with the exception of Uranus which has its rotational axis on the ecliptic plane indicating a 
possible collision with substantial mass loss. 

The fact that the boundary condition at z max for equilibrium with q(z max ) = differs 
from that for collapse with q(z max ) ~ 1 should not cause concern. Any self-similar solution, 
either Eulerian or Lagrangian, relies on writing the temporal and spatial dependences of any 
dynamic variable in separable form. This excludes the initial period of evolution of the real 
system, which is not amenable to similarity treatment. As a result, the self-similar initial 
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time t = is not the initial time t = of the real collapsing system. Thus, the homologous 
boundary condition at r = need not be the boundary condition of equilibrium. 

For q(0) ^> 1 and with a smaller a, the first minimum plunges to negative q as in Fig. 4, 
and becomes positive again at z Q = 4.5. As a result, the external part becomes disconnected 
physically from the central core, and a cavity is formed between the central core and the 
external envelope. The envelope will collapse under the gravitational field of the core plus 
the self-gravit y. Such a configuration is mo st relevant to pre-supernova ( magnet ar- in- a- 
cavity model (jUzdensky fc MacFadyen 1 120071 ) ). The envelope of this disconnected system is 
described by 

2 dq iN ( M* 



W^ = -( 7 -l),^ + M-l) , (24) 



M{z) = / q 3 (z')3z' 2 dz' 



3 1 

q 

where the integral M(z) starts at a disconnected z Q > 0. With 7 = 4/3 and neglecting the 
self-gravity M(z), we have the monotonic envelope profile of 

6a?W = 6a,(; ) - g - l) (I - i) . (25) 
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5. Model Parameters and Sonic Surfaces 

This model has two parameters M and a to normalize time t and radial coordinate 
77, that need to be determined. To obtain these two im portant parameter s, we begin with 



observations of a collapsed configuration, such as Swift f jSwift et al. II2005I ). Observing the 
central core dimension, we denote it by T\ which separates from the first minimum. This 
position is given by 

ri(r end ) = y{T en d)r)i = y(r end )ri(0) , (26) 

where T en d is the end time of the evolution function with y(r en d) vanishingly small, but 
nonzero. This y(T en d) has to be estimated independently. From the peak core density and 
the asymptotic envelope density, we can get the normalized g(0) from p(0) = g 3 (0). With g(0) 
determined by observations, we get the analytic solution of p(z) with a given a = C] /v'jf. 
We choose a = ot\ such that the inflection point z\ = arji = 1 to get 

1 = zi = arji = ari(0) = ^ n(r end ) . (27) 

With y(0) = 1 and y{r en d) estimated, we can construct backwards in time the position r\ at 
t = giving 

a = — . (28) 

m 

If the collapsed configuration allows the identification of the number of oscillations on the 
envelope, we can compare with the analytic p(z) to locate z max . By estimating the total 
mass M of the collapsed configuration, we then have 

M = M z 3 max , (29) 



which gives Mq as the mass in unit volume of z 3 = 1. With Mq and a determined, we go 

y s0 



back to the parameter a — a± to get C] = aiGM a. We now calculate 



M(z max ) = / p(z')3z' 2 dz' = z 3 max . (30) 



Since p(z) oscillates about unity, M(z) needs to be integrated in z. However, M{z max ) can 
be evaluated easily as above, since the collapse is merely a redistribution of densities of an 
initially uniform sphere with p(z) = 1, thus giving the second equility of the above equation. 
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In the Eulerian formulation, the singular sonic surfaces divide the (x,v) two-parameter 
space in subsonic and supersonic regions, where the flow is subsonic/supersonic relative to 
the similarity profile (IWhitworth fc Summers II1985I ). Complete collapse solutions have to 
cross these regions on the sonic lines. In our Lagrangian formulation, the fluid velocity, 
which is not a variable, and the sound speed are 



dy rnn/r a/2 d y 



a, = ^ 



p 



1/2 



JPo 
Po 



1/2 



dt 



dr 
C s qC s (z) 



(31) 
(32) 



We define the time dependent sonic function 

1 



S(z,t) = — 



dy 



{2a)V 2 C s {z) dr 



(33) 



With q(0) = 5 and a = 5, the spatial part S(z) is shown in Fig. 5 with oscillations in z. 
At r = with dy/dr = 0, the time dependent sonic function S(z, r) = everywhere. As 
t > and dy/dr increases, S(z,t) rises up and will cross the horizontal line of unit value 
at some point z s where S(z s ,r) = 1. This is the time dependent sonic surface z s (r). As r 
increases, this sonic surface moves relative to the similarity profile towards the center. Since 
dy/dr is related to the fluid velocity v = r/dy/dr, the movem ent of z x (r) amounts to the flui d 



velocity or flow speed with respect to the similarity profile (IWhitworth fc Summers Ill985l ). 



We, therefor e, have a picture of sonic surface which is compatible to that of the Eulerian 



formulation ( lLou fc Wang 1120061 ) . Nevertheless, in the Lagrangian formulation, v is directly 



incorporated in the evolution function, and the spatial structures are described by Eq. 
which is analytic everywhere in the one dimensional parameter space 77. There is no need of 
sonic lines to cross subsonic/supersonic regions. 
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6. Self-Organization and Conclusions 

Contrary to laboratory self-organized phenomena (lHasegawa l ll985h . astrophysical events 
given enough time with specific patt erns of outcome are sel f-organized, by defin ition. Grav i- 
tational collapses, extragalactic jets ( Tsui fc Serbeto 20071 ). planetary nebulae ( Tsui 2008 ). 
and supernovae are some of the examples. If, and only if, self-similar description is able to 
represent the astrophysical configurations in late times, it could only mean that these self- 
similar solutions are the attractors of dynamic time evolutions from a large set of random 
initial conditions. For the specific case of gravitational collapse, the traditional Eulerian 
homologous formulation, where the fluid velocity v is one of the variables, has the time de- 
pendence of each variable expressed in a given power of a linear time scaling. The spatial 
dependence is described by a set of nonlinear differential equations. The sonic surfaces and 
sonic lines are crucial in constructing complete solutions in the (x, v) two-parameter space. 

We have presented a Lagrangian homologous formulation, where the fluid velocity v is 
not a variable, but is derived from the evolution function y(t). The fact that v is not a variable 
is extremely important, because the parameter space is now reduced to one dimensional in rj. 
In particular, no sonic lines are needed to cross subsonic/supersonic regions. The crossing is 
embedded in the spatial structures. The time dependence of each variable is expressed in a 
given power of the evolution function, which is solved consistent to similarity. This evolutio n 
function y(t) agrees with the well established parametric solution of Mestel (IMestel Ill965l ). 
The spatial structure in the Lagrangian coordinate rj is described by an equation similar 
to that derived by Goldreich and Weber (jGoldreich fc Weber Ill980l ) with solutions covering 
the central core and the external envelope. In particular, the spatial solutions generate 
configurations that are relevant to planetary system in protostar and cavity in pre-supernova. 
The most important points of this Lagrangian description are that the evolution function 
is not limited to a linear time scaling, and the spatial solutions appear to be compatible to 
physical systems. These factors greatly enhance the scope of homologous treatment. We 
have the believe that the Lagrangian homologous description might be the mathematical 
tool to describe some self-organized astrophysical phenomena. 
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Fig. 1. — The evolution function y in thin line and its negative time derivative —dy/dr in 
thick line are plotted as a function of the normalized time r. 
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Fig. 2. — The equilibrium profile of density p = q 3 in thick line and q in thin line are shown 
as a funtion of the normalized distance z with q(0) = 2 and a = 100 giving z max = 40. 
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Fig. 3. — The collapsing profile of density p = q 3 in thick line and q in thin line are shown 
as a funtion of the normalized distance z with q(0) = 5 and a = 5. The peak density around 
z = is off scale. 
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Fig. 4. — The collapsing profile of density p = q 3 in thick line and q in thin line are shown 
as a funtion of the normalized distance z with q(0) = 20 and a = 1 giving zq = 4.5. The 
peak p and q around z = are off scale. 
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z 



Fig. 5. — The spatial part of the sonic function S(z) in thick line and the time dependent 
sonic function S(z, r) in thin line at the moment of dy/dr = 0.5 are plotted as a funtion of 
the normalized distance z with g(0) = 5 and a = 5. The sonic surface is at z s = 3.95. 



